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We present a new type of the Hamiltonian replica-exchange method, in which not temperatures 
, but the van der Waals radius parameter is exchanged. By decreasing the van der Waals radii 

that control spatial sizes of atoms, this Hamiltonian replica-exchange method overcomes the steric 
restrictions and energy barriers. Furthermore, the simulation based on this method escapes from 
the local-minimum free-energy states and realizes effective sampling in the conformational space. 
\ ■ We applied this method to an alanine dipeptide in aqueous solution and showed the effectiveness of 

■yf^ ' the method by comparing the results with those obtained from the conventional canonical method. 
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I. INTRODUCTION 



Effective samplings in the conformational space by Monte Carlo (MC) and molecular dynamics (MD) simulations 
are necessary to predict the native structures of proteins. In the conventional canonical-ensemble simulations [l|, 0, 
however, it is difficult to realize effective samplings in complex systems such as proteins. This is because 
the usual canonical-ensemble simulations tend to get trapped in a few of many local-minimum states. To overcome 
these difficulties, various generalized-ensemble algorithms have been proposed (for reviews, see, e.g., Refs. 0, [![)■ 

The replica-exchange method (REM) is one of the most well-known methods among the generalized-ensemble 
algorithms (see Ref. for the MD version) . It is easier to implement than the multicanonical algorithm [ll|, Q3 ; 
which is also one of the most well-known generalized-ensemble algorithms (see Refs. [13|, [l4j for the MD version), 
because we do not have to determine a p robability weig ht factor in advance in the REM. In the multicanonical and 
i similar algorithms [HI, [3 l2(j| . l2ll I22I I23I l24j | we employ non-Boltzmann weight factors as the probability 

weight factors. These non-Boltzmann weight factors are not a priori known and have to be determined by tedious 
procedures. On the other hand, the usual Boltzmann weight factor is employed in REM, and therefore it is not 
necessary to determine the non-Boltzmann weight factor. The REM uses non-interacting replicas of the target system 
with different temperatures and realizes a random walk in temperature space by exchanging the temperatures of pairs 
of replicas. Accordingly, the simulation can avoid getting trapped in local-minimum free-energy states. 

For large systems such as proteins in aqueous solution, however, the usual REM has a difficulty. We need to increase 
the number of replicas in proportion to 0(/2), where / is the number of degrees of freedom Q. Large biomolecular 
systems, therefore, require a large number of replicas in the REM and hence huge amount of computation time. In 
order to overcome this difficulty, it was pointed out that the number of required replicas can be greatly decreased if 
only the parameter exch ang es are performed in the multi-dimensional replica-exchange method (MREM) [2oT | without 
temperature exchanges |26| . In MREM replica exchanges in temperature and/or parameter in the potential energy 
are performed. MREM is also referred to as the Hamiltonian replica-exchange method [26j |. and in the present article 
we use the latter terminology. 

When we perform simulations of a protein in explicit water solvent, most of the degrees of freedom is occupied 
by water molecules. In order to predict the native structure of a protein, for instance, we would like to sample 
effectively the conformational space of the protein rather than water molecules. As an application of the Hamiltonian 
REM, therefore, Berne and co-workers performed simulations of the peptide in explicit water solvent, in which the 
scales of the only interactions related to the protein are varied [27| . They could achieve effective samplings in the 
conformational space of the peptide and saved CPU cost in comparison with the usual REM. Moreover, another 
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application of the Hamiltonian REM was reported by Kannan and Zacharias (28|. They focused on the backbone 
dihedral angles of the peptides and added biasing potential energy to the backbone dihedral angles. They could also 
achieve effective samplings in this space with the biasing potential energy. However, these biasing potential energy 
terms are complicated functions and highly dependent on force fields. 

In this article, we propose a new type of Hamiltonian REM where we exchange the scaling factor of the van der 
Waals radius of solute atoms in the interaction terms among solute atoms only. By reducing this scaling factor, the 
steric hindrance among solute atoms will be reduced and wide conformational space can be explored. We applied this 
method to the sytems of an alanine dipeptide with explicit water molecules and tested the effectiveness of the method 
by comparing the results with those from the conventional canonical method. 

In Section [H] we describe the new Hamiltonian REM. We give the details of the new Hamiltonian REM and 
canonical simulations that we performed in Section |TTTJ The results are presented in Section [IV] Section [V] is devoted 
to conclusions. 



II. METHODS 



A. Hamiltonian REM 

We first give a general formulation of the Hamiltonian REM [25j . 

We consider a system of N atoms with their coordinate vectors and momentum vectors denoted by q = {q 1 , • • • , q N } 
and p = {pi, ■ ■ ■ ,Pjv}i respectively. The Hamiltonian H\ in state x = (q,p) is given by the sum of the kinetic energy 
K and potential energy Ey. 

H x (x) = K(p)+E x (q) . (1) 

Here, we are explicitly writing (or introducing) a parameter of interest in the potential energy as A. In the canonical 
ensemble at temperature T, each state x is weighted by the Boltzmann factor: 

W B (x) = e~ fiH ^ , (2) 

where the inverse temperature (3 is defined by (3 = (&b is Boltzmann's constant). 

The generalized ensemble for the Hamiltonial REM consists of M non-interacting copies (or, replicas) of the original 
system in the canonical ensemble at M different parameter values X rn (to = 1, • • ■ , M). We arrange the replicas so 
that there is always exactly one replica at each A value. Then there is a one-to-one correspondence between replicas 
and parameter values; the label i (i = 1, ■ • • , M) for replicas is a permutation of the label to (to = 1, • • ■ , M) for A m , 
and vice versa: 

i = i(m) = /(to) , _ 
to = m(i) = / _1 (i) , 

where /(to) is a permutation function of to and / _1 («) is its inverse. 

Let X — (xi^\--- , i^ M " J = {e^x)'"' ,^m(lv/)| stand for a "state" in this generalized ensemble. Here, the 

superscript i and the subscript m in Xm label the replica and the parameter, respectively. The state X is specified 
by the M sets of coordinates q^ and momenta of N atoms in replica i at parameter A m : 



• (4) 



Because the replicas are non-interacting, the weight factor for the state X in this generalized ensemble is given by 
the product of Boltzmann factors for each replica i (or at each parameter A m ): 



M M 



WWmPO = flex? {-/3H Xm(i) (x® (i) ) } = II ex P > ( 5 ) 

i—l m— 1 

where i(m) and m(i) are the permutation functions in Eq. ([5]). 

We now consider exchanging a pair of replicas in the generalized ensemble. Suppose we exchange replicas i and j 
which are at parameter values X m and A„, respectively: 

X = • • , x\\ , • • • , x\^ , • • • \ ► X — ^ • • • , x\jj , • • • , x'Q , • • • \ . (6) 
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The transition probability for this replica exchange process is given by the usual Metropolis criterion: 

w(X -> X') = min (l, ^mS ) = ^ (1 ' 6Xp( " A)) ' (7) 



where we have [2f 



A = (3 



(E Xm ( g M) - J5 Am ( ff W)) - (^ An ( ff W) - £ A „ ( g W))] . (8) 



Here, and stand for coordinate vectors for replicas i and j, respectively, before the replica exchange. Note 
that we need to newly evaluate the potential energy for exchanged coordinates, E\ m (q^) and E\ n (q^), because E\ m 
and E\ n are in general different functions. We remark that the kinetic energy terms have canceled out each other in 
Eq. 

The Hamiltonian REM is realized by alternately performing the following two steps jH| : 

1. For each replica, a canonical MC or MD simulation at the corresponding parameter value A m is carried out 
simultaneously and independently for a certain steps with the corresponding Boltzmann factor of Eq. ^ for 
each replica. 

2. We exchange a pair of replicas i and j which are at the parameter values X m and A n , respectively. The transition 
probability for this replica exchange process is given by Eqs. ([7]) and ©. 

Finally, we remark that in order to further enhance sampling, we can always introduce a temperature-exchange 
process together with the above parameter exchange [25j . 



B. van der Waals Replica- Exchange Method 

We now describe our special realization of the Hamiltonian REM, which we refer to as the van der Waals Replica- 
Exchange Method (vWREM). 

We consider a system consisting of solute molecule(s) in explicit solvent. We can write the total potential energy 
as follows. 

E x (q) = E p {q p ) + E ps {q p , q s ) + E s (q s ) , (9) 

where E p is the potential energy for the atoms in the solute only, E ps is the interaction term between solute atoms 
and solvent atoms, and E s is the potential energy for the atoms of the solvent molecules only. Here, q = {q pi q s }, 
where q p and q s are the coordinate vectors of the solute atoms and the solvent atoms, respectively, and denoted by 

<7 p = jqfi, • ■ ■ , q Np | and q s = |q A r p+1 , ■ ■ • , q N ^. {N p is the total number of atoms in the solute.) 

We are more concerned with effective sampling of the conformational space of the solute itself than that of the 
solvent molecules. The steric hindrance of the solute conformations are governed by the van der Waals radii of each 
atom in the solute. Namely, when the van der Waals radii are large, the solute molecule is bulky and we have more 
steric hindrance among the solute atoms by the Lennard- Jones interactions, and when it is small, the solute molecule 
can move more freely. We thus introduce a parameter A that scales the van der Waals radius of each atom in the 
solute by 

<Tkt > Xa M , (10) 

and write the Lennard- Jones energy term within E p in Eq. ^ as follows: 

where r^i is the distance between atoms fc and I in the solute and eki and o~m are the corresponding Lennard- Jones 
parameters. The original potential energy is recovered when A = 1, and the steric hindrance of solute conformations 
is reduced when A < 1. Note that this is the only A-depcndent term in E\ in Eq. (J9j> . 

We prepare M values of A, A m (m = 1, • • • , M). Without loss of generality, we can assume that the parameter 
values are ordered as Ai < • • • < Xm- The vWREM is realized by alternately performing the following two steps: 
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1. For each replica, a canonical MC or MD simulation at the corresponding parameter value A m is carried out 
simultaneously and independently for a certain steps with the corresponding Boltzmann factor of Eq. (J5]) for 
each replica. 

2. We exchange a pair of replicas i and j which are at the neighboring parameter values A m and X m +i, respectively. 
The transition probability for this replica exchange process is given by Eq. ([7]), where A in Eq. ([8]) now reads 



A = /3 



(Vx m - V\ m ($)) - (V\ m+1 _ y w ($))] . (12) 



Here, V\ is the Lennard- Jones potential energy in in Eq. (|11|) among the solute atoms only. 

Note that because the A dependence of E\ exists only in V\, the rest of the terms have been cancelled out in Eq. {SJ. 

We see that Eq. ([T2")) includes only the coordinates g p of the atoms in the solute only and is independent of the 
coordinates q s of solvent molecules. Because N p <C N usually holds, the difficulty in the usual REM that the number 
of required replicas increases with the number of degrees of freedom is much alleviated in this formalism. 



C. Reweighting Techniques 

The results from Hamiltonian REM simulations with different parameter values can be analyzed by the reweighting 
techniques [2^ . |3C| . Suppose that we have carried out a Hamiltonian REM simulation at a constant temperature To 
with M replicas corresponding to M parameter values A m (m = 1, • • • , M). 

For appropriate reaction coordinates £i and £2, the canonical probability distribution Pt,\(£,i,(,2) with any param- 
eter value A at any temperature T can be calculated from 

M 

J2 (9m)- 1 N m (E Xl ,--- ,E XM ;^,^)e-^ 

Pr,x(tu&)= E — m > ( 13 ) 

E (dmV 1 rime 1 



M 



and 



-/t ,a„ 



E^o,A m (6,6)- (14) 

6,6 



Here, g m = 1 + 2r m , r m is the integrated autocorrelation time with the parameter value A m at temperature To, 
N m (E\ 1 , • ■ • , E\ M ] £1, ^2) is the histogram of the M -dimensional energy distributions at the parameter value A m and 
the reaction coordinate values (^1,^2), which was obtained by the Hamiltonian REM simulation, n m is the total 
number of samples obtained at the parameter value A m . Note that this probability distribution is not normalized. 
Equations (fl"3|) and (fl4|) are solved self-consistently by iteration. For biomolecular systems the quantity g m can safely 
be set to be a constant in the reweighting formulas [30] , and so we set g m = 1 throughout the analyses in the present 
work. Note also that these equations can be easily generalized to any reaction coordinates (£1,62, ■ ■ • )■ 

^From the probability distribution Pt,\(£i,&) ha Eq. (flU)) , the expectation value of a physical quantity A with any 
parameter value A at any temperature T is given by 

<^> T . A = - • (15) 

We can also calculate the free energy (or, the potential of mean force) as a function of the reaction coordinates £1 
and £ 2 with any parameter value A at any temperature T from 

F T ,\{&>&) = -fcBTlnP T ,A(6,6) ■ (16) 

By utilizing these equations, therefore, we can obtain various physical quantities from the Hamiltonian REM 
simulations with the original and non-original parameter values. We remark that although we wrote any T in 
Eqs. (fT3")l . (jTSJ) , and ((THJ) above, the valid value T is limited in the vicinity of T . We also need the T-exchange process 
in order to have accurate average quantities for a wide range of T values. 
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III. COMPUTATIONAL DETAILS 



In order to demonstrate the effectiveness of the present Hamiltonian REM, namely, vWREM, in which we exchange 
pairs of the van der Waals radius parameter values, we applied the vWREM MD algorithm, which we refer to as 
the vWREMD, to the system of an alanine dipeptide in explicit water solvent and compared the results with those 
obtained by the conventional canonical MD simulation. The N-terminus and the C-terminus were blocked by the acetyl 
group and the N-methyl group, respectively. The number of water molecules was 67. The force field that we adopted 
was the AMBER parm96 parameter set [3l| , and the model for the water molecules was the TIP3P rigid-body model 
[32| . The vWREMD and canonical MD simulations were carried out with the symplectic integrator with rigid- body 



conditions. The electrostatic potential energy was calculated by the Ewald method, and we employed the minimum 
image convention for the Lennard- Jones potential energy. The time step was taken to be 0.5 fs. 

In the vWREMD simulation, we needed only four replicas [M — 4). That is, we employed four different parameter 
values A m (m = 1, • • ■ ,4). In TableQ]we list the values of these parameters. The original potential energy corresponds 
to the scale factor A4 = 1.0, and in the canonical MD simulation we employed this original parameter value of 1.0, 
The temperature of the system To was set to be 300 K for all the replicas in the vWREMD and in the canonical MD 
simulation. Moreover, the initial conformations were also the same for all the simulations, and the initial backbone 
dihedral angles cf> and ip of the alanine dipeptide were set (<fi, ip) — (180°, 180°) as shown in Fig. [T] The total time 
of the MD simulations were 2.5 ns per replica for the vWREMD simulation and 2.5 ns for the canonical simulation, 
including equilibration for 0.1 ns. The trajectory data were stored every 50 fs. The replica exchange was tried every 
250 fs in the vWREMD simulation. 



IV. COMPARISONS OF THE VWREMD SIMULATION WITH THE CANONICAL MD SIMULATION 

We first examine whether the exchanges of pairs of the parameter values were realized sufficiently in our vWREMD 
simulation. In Table |TT] we list the acceptance ratios of replica exchange of the parameter values in the vWREMD 
simulation. These acceptance ratios are large enough (> 40 %) except for the pair of A4 and Ai. Because the difference 
of these parameter values is much larger than those of other pairs, this low acceptance ratio is expected, and it does 
not affect the REM performance. Figure shows the time series of the parameter set number m in A m which each 
replica visited. This figure shows that random walks in the parameter space were realized in the vWREMD simulation. 

Figures [3] and 0] show the time series of the backbone-dihedral angles tj> and ip for each replica. Figure [3] indicates 
that it is difficult to sample the range of the dihedral angle cp between 90° and 180°. This is because the steric 
restriction between the Cp atom of the alanine and the O atom of the acetyl group prevent the rotation in this range. 
For instance, Fig. [5] is the snapshot of the alanine dipeptide which corresponds to the structure at 2,469 ps in Replica 
3. In this structure the angle (p is 148.8°, and the distance between the atom of the alanine and the O atom of 
the acetyl group is 2.77 A. For such a short distance, the two atoms collide with each other. Therefore, the sampling 
among the range of 90° < cj) < 180° in the vWREMD simulation was rare, although the scale factor A m for the van 
der Waals radii was lessened. In order to sample this range frequently, it is necessary to employ a much smaller scale 
factor than the present case. However, conformations among this range have quite high potential energy due to the 
collisions between the Cp atom of the alanine and the O atom of the acetyl group, and it is not so important to 
sample this range at room temperature. On the other hand, the vWREMD simulation realized effective samplings 
with respect to ip in all the replicas as shown in Fig. [U This is because the van der Waals radius of the H atom of the 
N-methyl group that have the covalent bond with the N atom of the N-methyl group is small, and steric restrictions 
between the Cp atom of the alanine and the H atom are less than those between the atom of the alanine and the 
O atom of the acetyl group. 

Figures [5] and [7] show the time series of the backbone-dihedral angles <p and ip for each parameter value. For 
comparisons, we also show those from the conventional MD simulation with the original parameter value of A = 1. 
These figures show that the smaller the scale factor A m of the van der Waals radii is, the more efficient the sampling 
in the backbone-dihedral- angle space is. This is because the steric restrictions are reduced by lessening the scale 
factor and the energy barriers caused by the steric restrictions are decreased. Comparing Fig.^d) with Fig.[S]Je) and 
Fig- Eld) with Fig. [T^e), moreover, the vWREMD with the original parameter value A4 sampled the dihedral-angle 
space more effectively than the usual canonical MD simulation and did not get trapped in the local-minimum free- 
energy states. In other words, the canonical MD simulation with the original parameter value could not overcome 
energy barriers caused by the steric restrictions and got trapped in the local-minimum free-energy states. Therefore, 
effective samplings in the dihedral-angle space cannot be realized in the usual canonical MD simulations. 

The logarithm of the probability distributions, lnP((p,ip), with respect to dihedral angles <p and ip were obtained 
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from the vWREMD simulations. Figures [8] and [9] show lnP(<fi, ip) for each replica and those for each parameter value, 
respectively. From Fig.[5]we see that there are two regions in the vicinites of (</>, ip) = (0°, 0°) and (<p, ip) — (0°, 180°) in 
which the samplings were rare except the range of 90° < <p < 180°. This is because in the vicinity of ((/>, ip) — (0°, 0°), 
the O atom of the acetyl group and the H atom of N-methyl group collide with each other as in a structure shown in 
Fig. fTOT a). In the region of the neighborhood of (<p, ip) = (0°, 180°), the O atom of the acetyl group and the O atom 
of the alanine also collide as in a structure shown in Fig. rTOT b). It is obvious that steric restrictions in these regions 
were reduced and that the energy barriers were decreased by decreasing the scale factor as shown in Fig. [5] However, 
the effects of reducing energy barriers around (4>,ip) = (0°, 180°) are less than those around (<j>,ip) — (0°,0°). This 
is because the O atom of the acetyl group and the O atom of the alanine have negative charges, and repulsive forces 
caused by these charges act on both atoms. Therefore, by reducing the atomic charges as well as the van der Waals 
radii, more effective sampling may be realized in the conformational space. 

Figure QT] shows the free-energy landscapes at To = 300 K with respect to the backbone-dihedral angles </> and ip 
with the original parameter value (A = 1). The free-energy landscape in Fig. [TTT a) was calculated from Eq. (fT6|) 
by the reweighting techniques in Eqs. (fT3|) and (fl"4|) . The free-energy landscapes in Fig. fTITb) and Fig. [TlTc) were 
obtained from the raw histogram with the original parameter value A4 in the vWREMD simulation and the raw 
histogram in the canonical MD simulation, respectively. The free-energy landscape of the conventional canonical MD 
simulation is inaccurate due to insufficient sampling in the backbone-dihedral-angle space as shown in Figs. [6^e) and 
[T^e). The free-energy landscape obtained by the reweighting techniques in Fig. [TTT a) shows a better statistics even 
at free-energy barriers among the local-minimum states in comparison with that obtained from the raw histogram 
in Fig. [TTT b). This is because the information of all the other parameter values can be reflected by the reweighting 
techniques. Therefore, more accurate free-energy landscape can be obtained by the reweighting techniques. 

V. CONCLUSIONS 

In this article, we introduced a new type of Hamiltonian REM, the van der Waals REM (vWREM) in which the 
scale factor of the van der Waals radii of the solute atoms is exchanged only in the Lennard- Jones interactions among 
themselves. The steric hindrance due to the Lennard-Jones repulsions can be reduced by this method. Accordingly, 
the vWREM simulation can realize effective sampling in the backbone-dihedral-angle space without getting trapped 
in local-minimum free-energy states in comparison with the conventional canonical MD simulation. Employing the 
reweighting techniques, furthermore, we can obtain accurate free-energy landscape. 

Although we considered only exchanges of the scale factor of the van der Waals radii in this article, this idea can 
be extended to other parameters. For example, the scale factor of partial charges of solute atoms so that we can also 
realize even more efficient sampling in the conformational space. The generalization of the formalism in Sec. [IT] to 
other parameters is straightforward including the reweighting techniques. Moreover, these formalisms are independent 
of the degrees of freedom of solvent molecules. Therefore, these algorithms can be easily applied to large biomolccular 
systems. 
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FIG. 1: The common initial conformation of the alanine dipeptide. The figure was created with RasMol [39l |. 
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FIG. 2: Time series of the label m of the parameter \ m (m = 1,2,3,4) in (a) Replica 1, (b) Replica 2, (c) Replica 3, and (d) 
Replica 4 during the vWREMD simulation. 
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FIG. 3: Time series of the dihedral angles <j) m ( a ) Replica 1, (b) Replica 2, (c) Replica 3, and (d) Replica 4 during the 
vWREMD simulation. 
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FIG. 4: Time series of the dihedral angles ij) in (a) Replica 1, (b) Replica 2, (c) Replica 3, and (d) Replica 4 during the 
vWREMD simulation. 




FIG. 5: The snapshot of the alanine dipeptide at 2,469 ps in Replica 3 in the vWREMD simulation. The distance between 
the C/3 atom of the alanine and the O atom of the acetyl group is presented. The figure was created with RasMol [39l ]. 
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FIG. 6: Time series of the dihedral angle <j> with the parameter values (a) Ai, (b) A2, (c) A3, and (d) A4 during the vWREMD 
simulation, (e) Time series of the dihedral angle <j> during the conventional canonical MD simulation. 
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FIG. 7: Time series of the dihedral angle ip with the parameter values (a) Ai, (b) A2, (c) A 3 , and (d) A4 during the vWREMD 
simulation, (e) Time series of the dihedral angle ip during the conventional canonical MD simulation. 
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FIG. 8: The logarithm of the probability distribution \nP(<fi,ip) with respect to the dihedral angles cf> and ip for (a) Replica 1, 
(b) Replica 2, (c) Replica 3, and (d) Replica 4 in the vWREMD simulation. 



(a) 




-180 -90 90 180 

to 



(o) 




-180 -90 



90 180 



4>n 



(b) 




-180 



-180 -90 90 180 

4>C) 



(d) 




-180 



-180 -90 



90 180 



4>C) 



FIG. 9: The logarithm of the probability distribution \nP(4>, ip) with respect to the dihedral angles <f> and tp with parameter 
values (a) Ai, (b) A2, (c) A3, and (d) A4 in the vWREMD simulation. 
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FIG. 10: (a) The snapshot of the alanine dipeptide at 895.55 ps in Replica 1 in the vWREMD simulation. The structure has 
backbone-dihedral angles of (</>, ip) = (15.5°, 5.2°). The distance between the O atom of the acetyl group and the H atom of 
N-methyl group is given, (b) The snapshot of the alanine dipeptide at 2,465.65 ps in Replica 1 in the vWREMD simulation. 
The structure has backbone-dihedral angles of (4>,ip) = (—27.7°, 174.1°). The distance between the O atom of the acetyl group 
and the O atom of the alanine is given. The figures were created with RasMol [39l ]. 
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FIG. 11: The free-energy landscapes at To = 300 K with respect to the backbone-dihedral angles <f> and ip with the original 
parameter value A = 1. These were obtained from (a) the reweighting techniques applied to the results of the vWREMD 
simulation, (b) the raw histogram with the original parameter value A4 in the vWREMD simulation, and (c) the raw histogram 
in the canonical MD simulation. 



